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Clustering of Luminous Red Galaxies III: 
Detection of the Baryon Acoustic Peak in the 3-point Correlation Function 
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We present the 3-point function ^3 and Qa = C3/C2 for a spectroscopic sample of luminous red 
galaxies (LRG) from the Sloan Digital Sky Survey DR6 & DR7. We find a strong (S/N>6) detection 
of Qa on scales of 55-125 Mpc/h, with a well defined peak around 105 Mpc/h in all ^2, ^3 and Qa, 
in excellent agreement with the predicted shape and location of the imprint of the baryon acoustic 
oscillations (BAO). We use very large simulations (from a cubic box of L=7680 Mpc/h) to asses 
and test the significance of our measurement. Models without the BAO peak are ruled out by the 
Qi data with 99% confidence. This detection demonstrates the non-linear growth of structure by 
gravitational instability between z = 1000 and the present. Our measurements show the expected 
shape for Qz as a function of the triangular configuration. This provides a first direct measurement 
of the non- linear mode coupling coefficients of density and velocity fiuctuations which, on these large 
scales, are independent of cosmic time, the amplitude of fluctuations or cosmological parameters. 
The location of the BAO peak in the data indicates Q.m = 0.28 ± 0.05 and Q.b = 0.079 ± 0.025 
(for h = 0.70) after marginalization over spectral index (?is = 0.8 — 1.2) linear 61 and quadratic C2 
bias, which are found to be in the range: h\ = 1.7 — 2.2 and C2 = 0.75 — 3.55. The data allows a 
hierarchical contribution from primordial non-Gaussianities in the range Qs = 0.55 — 3.35. These 
constraints are independent and complementary to the ones that can be obtained using the 2-point 
function, which are presented in a separate paper. This is the first detection of the shape of Qs on 
BAO scales, but our errors are shot-noise dominated and the SDSS volume is still relatively small, 
so there is ample room for future improvement in this type of measurements. 
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I. INTRODUCTION 



The galaxy three-point function ^3 provides a crucial 
test and valuable statistical tool to investigate the ori- 
gin of structure formation and the relationship between 
galaxies and dark matter (see Bernardeau et al. [l| for 
a review). We will concentrate here on the reduced 3- 
point function Qs — ^3/^2 defined in Eql3]as the scaling 
expected from non-linear couplings (with Q^ ~ 1). Mea- 
surements of the three-point function and other higher- 
order statistics in galaxy catalogs have a rich history (eg 
Peebles and Groth Q, Fry and Peebles Q, Baumgart 
and Fry [1] Gaztanaga 0, Bouchet et al. @, Fry and 
Gaztanaga Q). In the past decade, three-point statistics 
have confirmed the basic picture of gravitational instabil- 
ity from Gaussian initial conditions (Frieman and Gaz- 
tanaga [a] , Jing and Boerner [9| , Frieman and Gaztanaga 
jl(j |. Feldman et al. [ll|). The connection between these 
observables and theoretical predictions is best done on 
large scales, where the physics (gravity) is best under- 
stood, but the surveys previously available were not large 
enough to have good statistics on sufficiently large scales. 

With the completion of large redshift surveys such as 
2dFGRS (CoUess et al. [13) and SDSS (York et al. [H) 
we expect measurement of higher-order statistics to pro- 
vide tighter constraints on cosmology (Colombi et al. 
[l4J . Szapudi et al. [1^, Matarrese et al. [lg|, Scocci- 
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marro et al. [13], Sefusatti and Scoccimarro [3). First 
measurements of the redshift space ^3 in the 2dFGRS 
(Gaztafiaga et al. [H) and SDSS (eg Nichol et al. ^) 
show good agreement with expectations (see also Kayo 
et al. [2l|, Nishimichi et al. [23, Kulkarni et al. [23|] and 
references therein). 

We will assume here that the initial conditions are 
Gaussian. Current models of structure formation pre- 
dict a small level of initial non-Gaussianities that we 
can neglect here. A popular parametrization of this ef- 
fect is to assume that initial curvature perturbations, 
given by the gravitational potential, $, are given by $ = 
'^l + Inl {^\- < *i >); where <^l is Gaussian and Jnl 
is a non-linear coupling parameter of order unity. This 
produces non-Gaussianities in the matter density pertur- 
bations at wavenumber k which, using the Poisson equa- 
tion, are suppressed by the square of the horizon scale 
kn = Ho/c, so that: Q3 - 3fNL{kH/kf-T{k)/D{a), 
where T{k) ~ 1 is the so called CDM transfer func- 
tion and D(a) is the growth factor. In our analysis 
{kn/kY — 10^'^ so that these type of primordial non- 
Gaussianities produce negligible contribution to Q^ for 
models with /jv^ ~ 1. A more detailed analysis of 
this will be presented elsewhere (see Sefusatti and Ko- 
matsu [2J| for a detailed forecast for this model). If 
the non-Gaussianities come from a non-linear coupling 
in the matter density field rather than in the gravita- 
tional potential, the resulting 3-point function will have 
a non-Gaussian contribution similar to that produced by 
non-linear bias C2 (ie see Eg IT^ below and Conclusions). 

The shape and amplitude of Q^ depends on galaxy 



bias, ic how galaxy light traces the dark matter (DM) 
distribution. This is both a problem and an opportunity. 
A problem because biasing can confuse our interpretation 
of the observations. An opportunity because one can try 
to measure the biasing parameters out of Q^. This idea 
was first proposed by Fry and Gaztanaga [25[ and has 
been applied to the 3-point function and bispcctrum of 
real data ( [1, M, M M, M 113, [28,] ) . 

This paper is the third on a series of papers on clus- 
tering of LRG. In the first two papers [23, l30| we stud- 
ied redshift space distortions on the 2-point correlation 
function. The reader is referred to these papers for more 
details on the LRG samples, simulations and the system- 
atic effects. Similar LRG samples from SDSS have al- 
ready been used by different groups to study the 2-point 
function (eg [3l|, [S^, [SJ, [3J, [3^) and found good agree- 
ment with predictions in the BAO scales, where density 
fluctuations arc at a level of few percent. This is encour- 
aging and indicates that this sample is large and accurate 
enough to investigate clustering on such large scales. 

In this paper we follow closely the methodology pre- 
sented in 3 previous analysis. Barriga and Gaztanaga 
[36| presented a comparison of the predictions for the 
two and three-point correlation functions of density fluc- 
tuations, ^2 and ^3, in gravitational perturbation the- 
ory (PT) against large Cold Dark Matter (CDM) sim- 
ulations. Here we use the same method and codes to 
estimate the clustering in simulations. Gaztafiaga and 
Scoccimarro J37| extend these results into the non-linear 
regime and focus on the effects of redshift distortions and 
the extraction of galaxy bias parameters in galaxy sur- 
veys. Gaztafiaga et al. [l^ apply this methodology to 
the 2dFGRS. Here we apply the very same techniques to 
the LRG data, so the reader is referred to these papers 
for more details. 

Kulkarni et al. [23| have also estimated <^3 using LRG 
galaxies from SDSS DR3, but focusing on smaller scales. 
We use DR6 which has 3 times the area (and volume) 
of DR3. We also use a volume limited sample and a dif- 
ferent estimator for the correlation functions and errors, 
focusing on the largest scales. 



II. THEORY 



to be a fair sample of possible realizations (see Peebles 
1980). A possible complication with this approach is the 
so call finite volume effects which result in potential es- 
timation and ratio biases (eg see [l|, [33|)- For our sam- 
ples we have checked using a large simulation that these 
potential estimation biases are small compared to the er- 
rors. This can be seen in Fig[2] below which shows a 
good agreement between predictions and mock simula- 
tions that have the same size than the SDSS sample that 
we used in our analysis. We use a simulation with about 
512 more volume than the data. We split this large sim- 
ulation into 512 subsamples and estimate the mean and 
error (from the variance) of the 2 and 3-point correla- 
tion functions in the 512 subsamples. We find that this 
mean agrees well, well within the error, with the cor- 
responding correlation estimated in the full simulation. 
This indicates that the finite volume effects arc negligible 
compared to errors (see also Fig. 5 in paper IV, |39|). 

It is convenient to define a Qa parameter as Groth and 
Peebles \M 



^3iri2,r23,ri3) 



(3) 



^3[ri2,r23,ri3} 

Cf = 6(^12)6(^23) +6(''12)6(»'13) + 6(''23)6(f'13), 

where we have introduced a definition for the "hier- 
archical" three-point function ^|^. Note that, by homo- 
geneity, the 3-point function 6 or Q3 can only depend 
on the distance between ri, r2 and r^. This involves 3 
variables that define the triangle formed by the 3 points. 
In principle Q3 could depend on the geometry and scale 
of the triangle. Here we will use two of the triangle sides 
''12, '"23 and /i, the co-sinus of the angle between ri2 and 
f23, which we call a. Thus Q3 = Q3{ri2,ri3, /J,). 

The Q3 parameter was thought to be roughly constant 
as a function of triangle shape and scale (Peebles [41| ) , a 
result that is usually referred to as the hierarchical scal- 
ing. Accurate measurements/predictions show that Q3 is 
not quite constant in any regime of clustering, although 
the variations of Q3 with scale and shape are small com- 
pared to the corresponding changes in 6 or 6 , specially 
at small scales. 



Definitions 



B. Qa and mode coupling 



The two and three-point correlation functions are de- 
fined, respectively, as 

6('"i2) = {S{n)6{r,)) (1) 

6(^12, r23,ri3) = {S{n)5{r2)S{r3)) (2) 

where S{r) = p{r)/p — 1 is the local density fluctuation 
about the mean p = (p), and the expectation value is 
taken over different realizations of the model or physical 
process. In practice, the expectation value is over differ- 
ent spatial regions in our Universe, which are assumed 



We will illustrate next how measurements of Q3 pro- 
vide a direct estimation of the non-linear mode coupling 
of density and velocity fluctuations. First consider the 
fully non-linear fluid equations that determine the grav- 
itational evolution of density fluctuations, (5, and the di- 
vergence of the velocity field, 9, in an expanding universe 
for a pressurelcss irrotational fluid. In Fourier space (see 
Eq. 37-38 in [l|): 



5 + e 



6 + He 



m, 



--Jdkidkia{ki,k2)9{ki)6{k2) (4) 

^H'^S ^ - J dkidkif3{ki,k2)e{ki)9{k2) 



where derivatives are over time dr — adt and H = 
dlna/dr. On the left hand side 6 = S{k) and 9 = 9{k) are 
functions of the Fourier wave vector k. The integrals are 
over vectors fci and fc2 constrained to k = fci2 = fc2 — fci. 
The right hand side of the equation include the non-linear 
terms which are quadratic in the field and contain the 
mode coupling fmictions: 



a = 



ki2 * ki 



(3- 



fcl2(fc2 *fcl) 
91.2 jL2 



(5) 



where "*" is the scalar product of the vectors (Eq.39 in 
[l|). This functions, a and /3, account for the mixing 
of Fourier modes. Note that this functions are adimen- 
sional and depend on the geometry of the triangle formed 
by the two wave vectors ki and k2 that contribute to 
k = k2 — ki. The scalar product indicates that mode cou- 
plings are larger when the modes are aligned. Physically 
this means that density and velocity gradients created by 
(non-linear) gravitational growth will tend to be parallel. 
Consider now the following perturbation expansions: 



^-E^» 



E 



(6) 



The first terms, 5i and 9i in the above series are the 
linear solution of the fluid equations Eq|31 where we ne- 
glect the quadratic terms in the equations. The linear 
term yields i5i = —9i for the first fluid equation in EqH] 
Combined with the second fluid equation yields the well 
known harmonic oscillator equation for the linear growth: 



^1 + HSi - -a„,HHi = 



(7) 



This equation is valid in Fourier or in configuration space. 
In linear theory each Fourier mode evolves independently 
of the other and they all grow linearly out of the initial 
fields with the same growth function, ie 5i = D{t)SQ, 
where Sq is the value of the field at some initial time and 
D(t) is the linear growth function, which is a solution 
to the above harmonic equation. If the initial field Sq is 
Gaussian then Si is also Gaussian. As shown above mode 
coupling is a non- linear effect. By construction, the next 
terms in the series of EqlSl ^2 and 02 , are assumed to be 
quadratic in the linear terms, ie 62 oc Sf. The solution for 
the second order can be found by just replacing the above 
expansion into the fluid equations keeping the second or- 
der terms in the equation. First order terms in the left 
side of the fluid equations in Eq[4] vanish by construction 
of the linear equation. We then find (see Eq.l56 in |l|): 



S2ik) 
F2 



dkidk2F2{ki, k2)6i{ki)di{k2) 



(8) 
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Thus the second order term in the expansion is S2 o< F2S1 

and contains the mode coupling information through F2 . 

Let's consider now the observables, which are the 2 and 

3-point correlation functions. Note that for an initially 



Gaussian field 61 is also Gaussian so that < i5^ >= 
and < Sf >= 3 < Sf >'^. At leading order the 2-point 
function ^2 is dominated by the linear evolution and the 
second to leading order term is zero because < S1S2 >= 
F2 < 5f >=0. Thus we have ^2(i) = £'(i)^6(0)- The 
linear contribution to ^3 is also zero and the leading order 
in ^3 comes from having 62 in one of the 3 points and Si 
in the other 2 points. For illustration, let us ignore for 
now the arguments of the 3 points, and see how things 
scale: 

6 -< S1S1S2 >oc F2 < 6t >= 3F2 < 6j >2oc F2e| (9) 

This is in agreement with the scaling ^3 = Q3C2 in. Eql3l 
We therefore have that Q3 at leading order is just propor- 
tional to F2 ■ This is exactly the case when we account for 
the arguments in the 3 points, as can be easily checked. 
The dependence of F2 on the triangular configuration, 
ie through EqlH is all contained in the corresponding 
triangular configurations of Q3, ie Q3 = (53(^12, ri3,/x). 
This shows the interest of measuring Q^: it provides a 
direct way to measure the mode coupling information in 
F2, which includes the fully non-linear coupling in EqlS) 
Higher order correlations just provide information about 
different combinations of coupling functions a and /?. 



C. Gravitational instability 

We can now address the following question: is large 
scale structure produced by gravitational growth from 
small Gaussian fiuctuations? We could answer this ques- 
tion by measuring ^2 (t) . If we knew some initial condi- 
tions ^2(0), we can then estimate D{t)'^ ~ £,2{t)/£,2{0) and 
compare to the linear solution of Eq[71 By measuring Q3 
we can directly compare to F2. This is independent of 
the linear test in D{t) and does not require knowledge of 
the initial conditions ^2(0). The result should in fact be 
independent of time. 

We can also test gravitational growth with indepen- 
dence of time or initial conditions by using the linear 
relation between density and velocities, i5i = —9i. But 
this is a test of linear evolution and requires measure- 
ments of the velocity field (this, in fact, is tested using 
redshift space distortions in Paper-I of this series). Q3 
only needs density fields and explores the non-linear sec- 
tor of gravity. 

If the origin of structure is non gravitational (or grav- 
ity is non-standard) but the initial conditions are Gaus- 
sian, then on dimensional grounds we would also expect 
£,2{t) = D{t)^£^2{0) and Q3 to be indepent of time. But 
both the amplitude and shape of Qt, as a function of tri- 
angle configuration could be quite different if the fiuid 
equations are different, either because of a non-standard 
cosmology or non-standard law of gravity (eg see |ll.[42|). 
In the standard case, the shape of F2 is mostly indepen- 
dent of cosmological parameters (eg ilm, J^a or fif,), time 
or the amplitude of fiuctuations [l|. As shown above, 
density and velocity gradients produced by non-linear 



evolution arc parallel which results in enhancement of 
clustering for elongated triangles. The exact amplitude 
as a function of triangular shape provides a finger print 
for non-linear gravitational growth. It is purely a non- 
linear effect that is not present in the initial conditions 
(which are assumed to be Gaussian with Q3 = 0). Un- 
fortunately, the prediction for Q^ depends not only on 
the mode coupling F2 but also on the slope of the ini- 
tial ^2(0). Models with relatively more large scale power 
(ie smaller slope) produce structures with larger coher- 
ence which give rise to more anisotropic structures and 
stronger shape dependence in Q3. Fortunately the shape 
of ^2 can also be estimated from data and one can then 
test the mode coupling predictions for F2 , as we will show 
below. 
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D. Shape dependence and BAO 

As mentioned above, there is a degeneracy in the shape 
of Qa between dark matter density fim and baryon den- 
sity ^B because they produce degenerate slopes in ^2- 
This degeneracy is broken by the presence of the BAO 
peak, but this requires a measurement at 100 Mpc/h 
scales. To illustrate this, consider a power law spectrum 
P{k) ~ Ak". In this case (see [l| and references therein): 



^3irl2,r23,^J.) 
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r23 
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(10) 



4(3-2(n-f-3) + (n + 3)V^) 

7 7l2 



ri2 
Un2)^{r23)+V 



where V stands for permutations of the indexes 123, and 
fj, is the co-sinus of the angle between ri2 and r23 1 which 
we call a. Here we will only show results as a function of 
a for fixed ri2 and r23. The above formula is only valid 
for a power-law spectrum. For CDM we will use a full 
calculation as explained below. 

Elongated or "collapsed configurations" are those with 
a ~ or a ~ 180 deg. We use the term "strong con- 
figuration dependence" when there is a significant differ- 
ence between the collapsed and the perpendicular config- 
urations. By "weak configuration dependence" we mean 
that Q3 is "hierarchical" (ie constant as a function of a). 
Q3{a) flattens with the decrease of the spectral index n 
and for larger n has a strong configuration dependence 
with a characteristic "V" shape. This corresponds to a 
larger probability of finding 3 points aligned, a direct 
consequence of gravitational infall that enhances the fila- 
ments that are characteristic of large scale structure both 
in simulations and real data. In a CDM spectrum this 
dependence on n translates into fiat Q3 (rounder struc- 
tures) at small scales progressively getting stronger V- 
shape as the effective n gets larger because of the CDM 
transfer function. At a fixed scale, the V-shape gets more 
pronounced when we increase flm or when we decrease 
r^B- This is just due to the effect of fim and fis in the 



FIG. 1: Perturbation theory predictions for the reduced 3- 
point function Q3 for different values of Qm and Qb , as labeled 
in the Figure. All cases are unbiased, except for dotted line 
which corresponds to Qm ~ 0.24 and SI;, = 0.03 with &i = 2 
and C2 = 1. This is for triangles with two sides fixed to 
ri2 = 33 ± 5 Mpc/h and ri3 = 88 ± 5 Mpc/h as a function of 
the interior angle a between this two sides. As a varies from 
0—180 degrees the third side of the triangle changes from 
r23 = 55 Mpc/h to r23 — 121 Mpc/h. For high values of fis a 
BAO peak (marked by the arrow) shows at a = 115 degrees, 
which corresponds to r2s = 106 Mpc/h. The precise location 
depends on the cosmology used. 



CDM transfer function, which changes the effective spec- 
tral index n. If we fix the shape of the P{k) spectrum, 
the values of flm and ^Ib have very little effect on Q3, as 
mentioned above. 

These points are illustrated in FigU] which shows per- 
turbation theory predictions for Q3. For this calculation 
we use Eq.(8) in Barriga and Gaztahaga [3y|, where dif- 
ferent scales and triangular configurations are also shown. 
Here we focus on the results of Q3 around the BAO scales. 
For high values of Qb a baryonic peak emerges in Q3 
and this peak could be used as a cosmological prove or 
to break the ^m — ^b degeneracy. The peak is present in 
all £,2, £,3 and Q3. As far as we know this is the first cal- 
culation illustrating how the BAO peak shows in Q3{a). 

Predictions in redshift space are harder to make (see 
Gaztaiiaga and Scoccimarro [37|), but on large scales sim- 
ulations show that redshift distortions in ^3 and ^| cancel 
out in Q3, and we do not expect deviations from PT tree- 
level predictions (see FigH] below). 



E. Biasing 

The value and shape of Q3{a) changes with galaxy 
bias. On very large scales we expect that galaxy fiuc- 
tuations Sq can be modeled as a local (but non-linear) 
function of the corresponding matter fiuctuations S, so 



that 5g — F[S]. For small fluctuations. (5 <C 1, wc can 
expand this local function as: 



fe - m -El -^ ' 



(11) 



where i = comes from the requirement that (5g) = 0. 
It then follows (see Fry and Gaztanaga [2^ , Frieman and 
Gaztanaga Q) that: 



^nr) 



?? 



T- (Q3 + C2) 



(12) 
(13) 



where C2 = &2/&1, and the ~ sign indicates that this is 
the leading order contribution in the expansion given by 
Eq. pip above. Thus, in general, the linear bias prescrip- 
tion is not accurate for higher-order moments even when 
(5 <C 1, the reason being that non-linearities generate non- 
Gaussianities of the same order as those of gravitational 
origin. The linear bias term bi can produce distortions 
in the shape of Qs, while the non- linear terms C2 only 
shifts the curve. This is illustrated by the dotted line 
in Fig[TJ It is therefore possible, but challenging, to use 
the shape of Qg in observations, when compared to the 
DM predictions, to separate &i from C2 in the above re- 
lation. This gives an estimate of the linear bias &i which 
is independent of the overall amplitude of clustering (eg. 
as). This approach has already been implemented for 
the skewness 5*3 [H, El, the bispectrum [I, [ll|, [H, ill 
or the angular 3-point function [lO[ and has been used to 
forecast future analysis |45j. 



III. SIMULATION AND ERRORS 

To check our codes and estimate errorbars, we have 
used a comoving output at z=0 of a MICE simulation, 
run in the super computer Mare Nostrum in Barcelona 
by MICE consortium (www.ice.cat/mice). The simula- 
tion contains 2048^ dark matter particles, in a cube of 
side 7680Mpc/h (which we call MICE7680), Qm = 0.25, 
Of, = 0.044, as = 0.8, n^ = 0.95 and h = 0.7. This simu- 
lation uses the power spectrum of Eisenstein and Hu [40| 
(EH from now on). To be consistent, the same approxi- 
mation has also been used for the predictions. We assume 
here that this EH fit is good enough approximation for 
the precision in our analysis, but we note that Sanchez 
et al. [43 found that on the BAO scale this approxima- 
tion is only accurate at the few percent level. We also 
use the EH fit with "no wiggles" which for each cosmo- 
logical model produces a match to the smoothed shape 
of the power spectrum without the BAO wiggles. This 
is equivalent to the removal of the BAO peak in the cor- 
relation function but keeping the same correlation as a 
based line. Note that we have chosen to displayed our 
results in terms of Qm and ^b for a fixed h = 0.7. 



The observed LRG galaxies are at a mean cosmic time 
of z = 0.3 rather than cosmic time z = in our simula- 
tions. In principle, this should result in a slightly larger 
amplitude for the clustering in the simulations, because 
it corresponds to a later cosmic time. But the simula- 
tions have a lower normalization, a^ = 0.80, than the 
value of as = 0.85 inferred for the LRG galaxies after 
accounting for the effect of bias (see Paper I) . This two 
effects partially compensate each other and the net result 
is that simulations have very similar clustering amplitude 
to that inferred in the LRG data. 

To simulate biasing we select groups of particles us- 
ing friend-of-friends with linking scale of 0.20. We find 
a total of 107 million groups with more than 5 particles 
(M > 1.87 X 10^^). These groups approximately cor- 
respond to DM halos when the number of par ticles in 
the group is larger than few tens of particles [48[ . When 
the number of particles is smaller than this, the group 
might not always correspond a virializcd DM halo. But 
in any case, these groups sample high density regions 
which are biased tracers of the dark matter distribution. 
We can produce different mock galaxy catalogs by choos- 
ing the group richness (ie mass or number of particles). 
Very massive groups are more biased and more rare than 
smaller groups. What we do here is to select the group 
mass cut-off to reproduce the galaxy clustering in the 
LRG observations. It turns out that when we do this 
we get a number density of groups which is quite simi- 
lar to the number density of LRG galaxies. Our mock 
catalogs are not realistic in the sense that we have not 
simulated the physics of galaxy formation. But wc only 
use these catalogs for error estimation. Errors depend 
on the statistical properties of the simulation and not on 
process that produce such statistics. For our purposes 
here, all we need from mock simulations is that they have 
approximately the same volume, bias b, number density 
and similar values of ^2 and ^3 that observations. Errors 
only depend on this quantities. 

We have also used a MICE simulation with 2048^ dark 
matter particles, in a cube of side 3072Mpc/h (which we 
caff MICE3072, same parameters as MICE7680) which 
has 15 times better mass resolution to check for mass 
resolution effects. We find very similar results in both 
cases, but obviously MICE7680 provides more sampling 
volume to estimate reliable errors. 

When we select groups with M > 2 — 4 x 10^^, both 
the clustering amplitude (61 ~ 1.9 — 2.2 for as = 0.8) 
and the number density (fi ~ 4 — 6 x 10^'"') are similar 
to the real LRG galaxies in our SDSS sample (the range 
reflects the fact that the actual number depend on LRG 
sample used). As the mock simulations are similar to 
the real data we will use the variance of clustering in the 
simulations to estimate errorbars for our analysis. The 
resulting errors from simulations are typically in good 
agreement with JK errors from the real data (see paper 
I for details). We will also check if we can recover the 
theory predictions for Qt, by using mock simulations with 
similar size as the real data. 





FIG. 2: Measurements of Qj, as in Fig[l] in dark matter (left panel) and groups (right panel) from the MICE7680 mock 
simulations. Mean and errors are estimated from 512 subsamples of l{Gpc/Kf' . This is compared to the PT predictions in 
real space (dashed line) and corresponding biasing predictions for groups (continuous lines) Triangles and squares correspond 
to measurements in real and redshift space respectively. Two dotted lines shows the result for two single llGpc/hf' realizations 
(#20 and #21), in redshift space. 



We have divided the big MICE7680 cube in 8^ sub- 
cubes, each of side 961Mpc/h and apply the LRG SDSS 
mask to obtain Nm = 512 LRG mocks from both dark 
matter and groups. From the Nm mock catalogs, we can 
estimate what we call the Monte Carlo (MC) covariance 
matrix: 



1 ^'" 



(14) 



fe=i 



where C(*)'^ i^ the measure of ^2, Ca or Qa iii the k-th 
mock simulation (A: = 1, ...Nm) and ^(i) is the mean over 
Nm realizations. The case i=j gives the diagonal error 
(variance). 

Fig (21 compares the mean and errors in the MICE7680 
mocks with the PT predictions (dashed lines) in real 
space. At these large scales, results in real space are 
almost identical to redshift space for both dark matter 
(left panel) and groups with M > 1.87 x 10^^ (right 
panel). The PT predictions and the biasing predictions 
work remarkably well. For bias we have used &i = 1.9 
and C2 = 0.2 (continuous hues). The value 61 = 1.9 is 
estimated empirically from the ratio of the 2-point func- 
tion of the groups to that of dark matter. This ratio is 
fairly constant for r > llMpc/h. The non- linear bias 
C2 produces a global shift that we have just fitted to the 
data. These estimated values of bi and C2 agree very 
well with halo model prediction [43 for the mass of these 
groups. We note that we show the mean of the 512 mocks 
with errors from the dispersion in redshift space, which 
are slightly smaller (10-30%) than in real space. Dot- 
ted lines show results from two of the 512 realizations, 
illustrating the strong covariance in (53(a). 
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FIG. 3: Here we show the first 6 eigenvectors in the SVD 
using 512 subsamples in MICEL7680. Bottom panel shows 
the first 3 principal components (dashed, continuous and long- 
dashed lines) and the top panel shows the next 3 components, 
ranked by amplitude. While the fim — Q.b constraints are 
degenerate using the first 3 components, this degeneracy is 
broken by using more eigenvectors because they can separate 
BAO features. 



We apply a x^ method to fit the models by inverting 
the covariance matrix C^ , as explained in Gaztafiaga and 
Scoccimarro [331. Before inverting Cy in Eq.(8), notice 
that the values of Cy are estimated in practice to within 
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FIG. 4: Contours of constant Ax = 1,2.3,4,6.2. and 9 obtained from fitting models to measurements in dark matter (left 
panel) and groups (right panel) in MICE7680 mock simulations using SVD with MC covariance matrix corresponding to a 
Survey of about 1 Gpc^/h^. The crossing lines show the input values in the simulations, while the cross shows the best fit 
value. At 1-sigma, best fit is in excellent agreement with the input. But at 2-sigma there is a secondary peak and a strong 
degeneracy in the Qb — ^m plane. 



a limited resolution in Eq[T31 Adj ~ \/w~- There- 
fore if the number of mocks TV^ is small or if there are 
degeneracies in dj , the inversion will be affected by nu- 
merical instabilities. In order to eliminate this problem, 
we perform a Singular Value Decomposition (SVD) of 
the matrix. By doing the SVD decomposition, we can 
choose the number of modes we wish to include in our 
X^ by effectively setting the corresponding inverses of the 
small singular values to zero. In practice, we work only 
with the subspace of "dominant modes" which satisfy 
A^ > •y/2/iVm which is the resolution to which we can es- 
timate the covariance matrix elements. Typically results 
converge after using a few singular values. Most of the 
times there is no gain in using more components which 
indicates that the effective number of degrees of freedom 
is smaller than the number of bins. Sometimes the fit 
gets bad when using larger number of components, which 
indicates instabilities in the covariance inversion. The re- 
sults presented here are quite robust. We find that they 
are basically the same when we use different covariance 
matrices, corresponding to samples with different ampli- 
tude of clustering or different shot-noise, within the range 
that roughly matches the data. In particular, results are 
quite similar when we use DM particles mocks instead of 
groups mocks to estimate Cy- . The absolute errors in Q^ 
are roughly the same for dark matter and groups, despite 
the difference in the amplitude. This give us confidence 
that what we find in real data is not an artifact of our 
analysis. 

Fig [3] shows the first 6 eigenvectors or principal com- 
ponents, which resemble an harmonic decomposition of 
Qaia), as proposed by Szapudi [50[. The first component 
is roughly a constant, like a monopole. the second compo- 



nent corresponds to the "V" shape mentioned above, like 
a quadrupolc, and the third (short-dashed line) is sim- 
ilar to a dipole. The next components resemble higher 
multipoles and can be combined to define localized struc- 
ture in Qsia). With the first two components it is not 
possible to break the Qb — ^m degeneracy, given the er- 
rorbars. We have checked this in both the LRG data and 
the simulations. As we increase the number of multipoles 
we can see how the degeneracy in the fls — ^m plane be- 
gins to break down. The BAO feature only shows in the 
higher components and breaks the Qb — ^m degeneracy 
for large values of Qb > 0.03. Lower values show no sig- 
nificant BAO peak, given the errors, and this results in 
a strong degeneracy in Qb — Qm ■ We have played with 
the models and find that even if we artificially reduce 
the errors we need more that 2 eigenvalues to have a Qb 
detection. 

In Fig|4] we show how well we can recover the values 
of Qm — Qb from the mocks shown in Figl^l We use 
9 singular values, as in the LRG data, but results are 
quite similar from 5 to 9 singular values. In the case 
of dark matter (left panel in Fig|4]), we do not include 
biasing parameters. In the case of groups (right panel) 
the results are marginalized for bi ~ 1.7 — 2.2 and C2 = 
0.0 — 5.0 . The best fit value is in excellent agreement 
with the input values both for Q„i — i^b and for bi — 
02- This illustrates that fitting bi and C2 from Q3 in 
groups could be used to roughly estimate the mass of the 
groups or galaxy clusters. But even when we recover well 
the input value, note how errors are large and produce 
degenerate values at 2-sigma level for the size of our mock 
sample. As we will see below, this is because Qb is low 
in the simulation. The volume of data in a single mock. 



of about 1 (Gpc/h)'^, is not large enough to break the 
flrn ~ ^b degeneracy for low values of fli,. This could 
be improved by using more shape configurations in Qa 
or smaller smoothing scales (we use cubic pixels with 11 
Mpc/h on a side). The second (weaker) minimum x^ that 
shows in FigHlfor larger values of fib — 0.09 corresponds 
the case where the BAO peak moves to smaller values 
of a and overlaps with the valley of the "V" shape in 
(53(a). The Qs amplitude is low and the relative error 
is large in this region. This allows a BAO peak to be 
compatible with the simulations. This will also affect 
our LRG measurements in observations, although in this 
case the relative error is much smaller thanks to the non- 
linear bias C2 which increases the overall amplitude of 
Q3{a) and allows a more significant detection of Qb- 



IV. DATA AND ANALYSIS 
A. Data Sample 

The luminous red galaxies (LRGs) are selected by color 
and magnitude to obtain intrinsically red galaxies in 
Sloan Digital Sky Survey (SDSS). See Eisenstein ct al. 
[5l| or http://www.sdss.org for a complete description of 
the color cuts. These galaxies trace a big volume, around 
lGpc^h~^, which make them perfect to study large scale 
clustering. LRGs are red old elliptical galaxies, which 
are usually passive galaxies, with relatively low star for- 
mation rate. Since they reside in the centers of big halos 
they are highly bias with 6 ~ 2, between regular galaxies 
and clusters. 

LRG's are targeted in the photometric catalog, via cuts 
in the (g-r, r-i, r) color-color-magnitude cube. Note that 
all colors are measured using model magnitudes, and all 
quantities are corrected for Galactic extinction following 
[52J . The galaxy model colors are rotated first to a basis 
that is aligned with the galaxy locus in the (g-r, r-i) plane 
according to: 

c^= (r-i) - (g-r)/4 - 0.18 

C||=0.7(g-r) + 1.2[(r-i) - 0.18] 

Because the 4000 Angstrom break moves from the g 
band to the r band at a redshift z ~ 0.4, two separate 
sets of selection criteria are needed to target LRGs below 
and above that redshift: 

Cut I for z < 0.4 

rPetro < 13.1 -H c,, / 0.3 

rPetro < 19.2 

\c^\ < 0.2 

mu^o < 24.2 mag arcsec~^ 

rpSF - Tmodel > 0.3 

Cut II for z > 0.4 

rPetro < 19.5 

\c^\ >0.45- (g-r)/6 
g-r > 1.30 + 0.25(r-i) 



iTiUbo < 24.2 mag arcsec ^ 

rpSF - r model > 0.5 

Cut I selection results in an approximately volume- 
limited LRG sample to z=0.38, with additional galaxies 
to z ~ 0.45. Cut II selection adds yet more luminous 
red galaxies to z ~ 0.55. The two cuts together result in 
about 12 LRG targets per deg^ that are not already in 
the main galaxy sample (about 10 in Cut I, 2 in Cut II). 
The radial distribution and magnitude-redshift diagrams 
for these galaxies are shown in Paper I |29l |. 

We k-correct the r magnitude using the Blanton pro- 
gram 'kcorrect' [53|. We need to k-correct the magni- 
tudes in order to obtain the absolute magnitudes and 
eliminate the brightest and dimmest galaxies. We have 
seen that the previous cuts limit the intrinsic luminosity 
to a range —23.2 < Mr < —21.2, and we only eliminate 
from the catalog some few galaxies that lay out of the 
limits. Once we have eliminated these extreme galax- 
ies, we still do not have a volume limited sample at high 
redshift. For the 2-point function analysis we account 
for this using a random catalog with identical selection 
function but 20 times denser (to avoid shot-noise) . The 
same is done in simulations. Computationally this is very 
time consuming (specially as there are 512 simulations) 
because it involves N'^ operations where A^ ~ 10^. For 
the 3-point, using the random catalogs would involve N^ 
operations which begins to be very challenging with cur- 
rent computer power. We will therefore use a different 
estimator based on a pixelization of the sample. This 
estimator is ideally match to a volume limited sample, 
where the full volume is equally sampled and there is no 
radial selection function (we still have angular mask and 
radial boundaries). 

For the 3-point function analysis presented in the paper 
we select a volume limited sample with —22.5 < Mr < 
—21.5 and z = 0.15 — 0.38 from the spectroscopic sam- 
ple of LRG in the SDSS DR6. We choose this particular 
sample because it is the best compromise between vol- 
ume and number density. There are about 40, 000 LRG 
galaxies in this sample (n ~ 4 x 10^^). 



B. Correlation functions 



We estimate the correlation functions with a fast al- 
gorithm described in some detail in Gaztahaga et al. 
[19J . Barriga and Gaztanaga [301 . This algorithm allows 
a fast calculation of two and three-point function for mil- 
lions of points. The first step is to discretizc the simula- 
tion box into Lsize^ cubic cells. We assign each particle 
to a node of this new latticed box using the nearest grid 
point particle assignment. We precalculate the list of 
relative neighbors to any given node in the lattice. To 
compute now the two-point and three-point correlation 



functions wc use: 



6 (''12) 
£.3{ri2,r23,ri3) 






E 



i,j,k 



1 



(15) 
(16) 



where i extends over all nodes in the lattice, j over the 
list of precalculated neighbors that are at a distance ri2 ± 
dr /2 from i, and k is over the neighbors at distance r23 ± 
dr/2 from j and ris ± (ir/2 from i. We take dr to be 
equal to the pixel size. 

There are two sources of errors in this estimation: a) 
Shot-noise which scales as one over the square root of 
the number of pairs or triplets in each bin b) sampling 
variance which scales with the amplitude of the correla- 
tions. It is easy to check that for the size and density of 
our sample the shot-noise term dominates over the sam- 
pling variance error. This has been checked in detailed 
by using DM simulations which have a large density and 
can be diluted to explore the shot-noise contribution to 
the error budget (see also Paper I). 

We choose cubical pixels of dr = 11 Mpc/h on the 
side. This is an adequate compromise to measure BAO. 
We want this number to be as large as possible to re- 
duce shot-noise and to average over many triangles in a 
fast way. On the other hand, wc need a good resolution 
to avoid loosing too much shape and BAO information. 
Once the pixel size is fixed to dr = 11 Mpc/h we are 
forced to use ri2 > 3dr which is the smaller distance 
that does not distort the Qa shape [3Q|. To reach the 
BAO scale we need ri3 > Sdr. So this fixes our choice 
of triangles. The 2-point function for this sample and 
pixelization is shown in Fig.Bll of Paper I [29J . 



C. Results 

After the release of DR6, Swanson et al. [53| provided 
mask information in a readily usable form, translating 
the original mask files extracted from the NYU Value- 
Added Galaxy Catalog [III , from MANGLE into Healpix 
format [S^ • Cabre and Gaztahaga [231 describe how they 
constructed a survey "mask" for LRGs and tested the im- 
pact of the mask on clustering measurements using mock 
catalogs. Using the same techniques, we have also ex- 
amined the correlation function of LRGs in DR7, which 
has become available since the submission of Cabre and 
Gaztaiiaga [20|. In Fig. [5] we plot a summary of the 
possible systematics in the estimation of the Q3 measure- 
ments. Our main result, that will be used for comparison 
with models, is shown as a shaded region, corresponding 
to the 1-sigma region. Results using different masks or 
the DR7 release are quite similar despite variations of 
^ 17% in the fraction of galaxies or area used in the 
different cuts. We also show comparison to the results 
using all galaxies (open squares), rather than just a vol- 
ume limited subsample. In this case we define density 
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FIG. 5: Estimations of Q3 as in Fig[l]froni observations using 
different samples and masks. Closed (red) squares and (red) 
shaded region show the main measurements and errors used 
in this paper, ie for a volume limited sample from DR6 with 
a magnitude range -22.5 < M < —21.5 and 0.15 < z < 0.38. 
Closed (blue) circles show the corresponding result in the DR7 
sample. Open and closed triangles show the DR6 results using 
the MANGLE mask of Swanson et al. (2008) with greater 
than 0.0 or 0.8 completeness fractions. Open squares use all 
magnitudes. The (black) continuous line correspond to our 
best fit model in Fig|6] 



fluctuations using the local mean density provided by the 
random catalogs, rather than the overall mean density as 
we do for volume limited catalogs. The result is noiser 
(because of the shot-noise in the random catalogs) but 
in very good agreement with the other measurements. 
Evolutionary effects (that change the mean density) or 
magnitude effects do not seem important in this sample. 
We conclude that the results are very robust to the sys- 
tematic variations that we have tried, and choose to use 
our default DR6 volume limited sample because this is 
the one that have been more tested and is the bases for 
our error analysis. 

Fig [6] shows a comparison of Q3(a) measurements in 
the LRG sample with models. There is a good resem- 
blance with what is expected from theory, ie compared 
to FiglU with a peak at a ~ 100 deg. which resembles 
much the BAO peak predicted by models. As our sig- 
nal is shot-noise dominated one may wonder if this peak 
could be produce by noise fluctuations. 

An indication that this signal is real comes from Fig|71 
Detection of the BAO peak using the hierarchical prod- 
uct ^1^ of 2-point function, as shown in Fig [71 is in excel- 
lent agreeincnt with previous detections [3l| and results 
over the very same sample in Paper I of this series [231 ■ 
The peak seems to be detected in all ^3, Q3 and ^2 and 
products. 

Left panel of FigO shows the x^ fit foi' ^b ~ ^m plane 
using 9 singular values with a total signal-to-noisc of 6.25. 
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FIG. 6: Values of Q3 as Fig[T] this time comparing three 
of the biased models (lines) to observations in LRG galaxies 
(symbols with errorbars). The total signal-to-noise in this de- 
tection is 6.25. Models have ilm = 0.26 and h = 0.7. with 
ilb = 0.03 (short-dashed lines) or £7;, = 0.06 (continuous line). 
Long-dashed uses EH fit with the no- wiggles (ie no BAO peak, 
but otherwise equal baseline spectrum) while the other lines 
includes the BAO peak in the models. We have marginal- 
ized over biasing parameters and spectral index and show the 
best fit in each case. The high f2(, = 0.06 model has a mini- 
mum ')^ —% (with 3 degrees of freedom) which is significantly 
smaller than the mini-mun x^ ~ 17 for the $!(, = 0.03 model. 
The BAO peak shows at a = 100 in both the model with 
Q,b ~ 0.06 and the data. The best fit model without BAO 
peak has x^ > Hi ie a probability smaller than 1% of being 
correct. 
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FIG. 7: Separate measurements of ^3 (top panel) and hierar- 
chical ^3 = ^2{ri2)^2ir23) + ^2{r2i)^2{ri3) + +^2iria)^2ir32) 
(bottom panel). The models are as in Fig[Sl ie Sl(, = 0.03 
(short-dashed line) and fib = 0.06 with (continuous line) and 
without wiggles (short-dashed lines), all with Qm ~ 0.26. In 
this case the prediction depends not only on the biasing pa- 
rameters, but also on the erg normalization. As can be seen 
in this figure, the model with large fls show a different shape 
and a BAO feature both in ^3 and ^2. Data follows the BAO 
predictions in both quantities, as well as in Q3 which is quite 
reassuring. 



We fix /i = 0.7 and marginalized over spectral index n^ = 
0.8 — 1.2 and biasing parameters 61 = 1.7 — 2.2 and C2 = 
0.0 — 5.0. In the right panel we show the same fit for the 
EH models with no- wiggles (ie no BAO peak). Note how 
the values of CIb — ^m become degenerate, as expected 
from our previous argument that the BAO peak helps 
to break the fls — ^m degeneracy. The best fit value 
without BAO peak is x^ = H as opposed to x^ = 6 
with the BAO peak. Thus, in relative terms the models 
with and without a BAO peak are between 2 and 3-sigma 
away. But note that in absolute terms, models without 
the BAO peak are ruled out with > 99% confidence. 

The WMAP5 best fit value (marked by a cross) is 
outside the (2D) 1-sigma join region, but inside the 
(2D) 2-sigma contours. The best fit value is for fi^ = 
0.079 ± 0.025 and we find that fls > 0.035 at 2-sigma 
level for any value of f2„i- 

If we fix the Qb — ^m to its best fit value, we find 
61 = 1.7-2.2 and C2 = 0.75-3.55. The value of the linear 
bias bl is in excellent agreement with what we found in 
Paper I by fitting redshift space distortions in the 2-point 
function, but the error here is larger. The value of the 
non-linear bias C2 — 2 is higher than the one we found in 



previous section for halos C2 — 0.2. This is not surprising 
as it is well known that more than one LRG can occupy 
a single halo, in which case C2 tends to be larger for a 
given bl [43| . Also note that a larger value of C2 makes 
the Q3 signal-to-noise larger in the LRG data than in the 
MICE7680 group mocks. This helps defeating the shot- 
noise and improves the significance of the BAO detection. 



V. CONCLUSIONS 

We have studied the large scale 3-point correlation 
function for luminous red galaxies from SDSS, and partic- 
ularly the reduced Q3 = Cs/CIj which measures the scal- 
ing expected from non-linear couplings. We find a well- 
detected peak at 105Mpc/h separation that is in agree- 
ment with the predicted position of BAO peak. This 
detection is significant since it is also imprinted in ^2 and 
^3 separately. We focus our interpretation in Q3 because 
it is a measure independent of time, as or growth fac- 
tor. It only depends on the shape of the initial 2-point 
function and the non-linear coupling of the gravitational 
interaction. Our result for Qa is in excellent agreement 
with predictions from Gaussian initial conditions. When 
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FIG. 8: Contours of constant Ax = 1,2.3,4,6.2. and 9 obtained from fitting models to data in Fig [6] using a SVD with 
covariance matrix from the MICE group simulations. Contours are marginalized over h\ — 1.7 — 2.2, Ua = 0.8 — 1.2 and 
C2 = 0.0 — 5.0 (we use h = 0.7). The crossing lines show the best WMAP5 fit. Right panels uses EH fit with the no- 
wiggles (ie no BAO peak, but otherwise equal shape) while left panel includes the BAO peak in the models. Best BAO fit is 
fim = 0.28 ± 0.05 and Ob — 0.079 ± 0.025 with x^ = 6 for 3 degrees of freedom (9 singular values minus 5 parameters in the 
fit). Probability for models with no BAO peak is less that 1% (x^ > H)- 



we use the Q3 data alone (with no fit to the 2-point func- 
tion) we are able to break the strong degeneracy between 
fim and flB (see Fig 15]). Our detection shows a clear pref- 
erence for a high value of Ub = 0.079 ± 0.025. This value 
is larger, but still consistent at 2-a with recent results of 
WMAP {Qb = 0.045). At 3-sigma level, the n^-^B be- 
comes degenerate. Models with no BAO peak are ruled 
out at 99% confidence level. 

We have used very large realistic mock simulations to 
study the errors. These simulations show that Qa is not 
significantly modified in redshift space so we can use real- 
space perturbation theory (see Fig[2]). This agreement 
also indicates that loop corrections are small on BAO 
scales [56| |. This analysis is independent from 2-point 
statistics, which tests the linear growth of gravity, since 
3-point statistics test the non-linear growth. A high value 
for f2s is also consistent with the analysis of the peak in 
the 2-point correlation function shown in Eisenstein et al. 
[sH and in Paper IV (Gaztanaga et al. [3^) of this series, 
which detect a slightly higher peak than expected. We 
have done all the analysis with just one set of triangle 
configurations, with fixed sides of ri2 = 33 ± 5.5 Mpc/h 
and ri3 = 88±5.5 Mpc/h, to center our attention to BAO 
scale. This is about optimal, but we notice that there is 
much more to learn from Qa, which will be presented in 
future analysis. Results on smaller scales are consistent 
with what we find here. 

Data is in excellent agreement with Gaussian initial 
conditions, for which Q3 ~ 0. But note that our 
quadratic bias detection C2 = 0.75 — 3.55 is degenerate 
with a primordial non-Gaussian (hierarchical) contribu- 
tion. Indeed the mean value of C2 — 2 seems larger in 



observations than in halo simulations, for which we find 
C2 — 0.2. We believe that this indicates that halos are 
sometimes occupied by more than one galaxy, which in- 
creases the effective value of C2 [43 • But if we are con- 
servative we can not rule out a primordial non-Gaussian 
contribution in the range Q^{Primordial) = C2 — C2 = 
[0.55,3.35]. 

Note that we have pixelized our data in cubical cells 
of side dr = \1 Mpc/h. This results in some lost of 
small scale information but allows for a very fast method 
to estimate 3-point function [36|. This is important for 
data, but more for simulations. In the MICE7680 sim- 
ulation there are close to iV = 10^^ particles. A brute 
force method to estimate 3-point correlation would re- 
quire N^ = 10'^^ operations, while our method based on 
pixels just needed 5 x 10^^ operations. 

Future surveys will be able to improve much upon our 
measurement here. A photometric survey with Az < 
0.003(1 -I- z) precision (corresponding to dr < 9 Mpc/h 
at z=0), such as in the PAU Survey [53| should have 
enough spatial resolution to measure Q'i{a) as presented 
in this paper (recall that we are binning our radial dis- 
tances in dr = 11 Mpc/h). Such survey could sample 
over 10 times the SDSS DR6 volume (ie to z=:0.9) with 
20 times better LRG number density (ie for L > L*). 
Fig|5] shows the forecast for such a survey, which we have 
simulated with the MICE7680 mocks in redshift space 
with a photo-z of Az < 0.003(1 -I- z) and for the same 
triangles as shown in our SDSS analysis. This is just il- 
lustrative, as we have not marginalized over biasing and 
other cosmological uncertainties. But note that the im- 
provenrent is substantial and shows the potentiality of 
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this method to constrain cosmological parameters and 
models of structure formation. 



FIG. 9: Same as Fig|4]for a future photometric Survey with 
photo-z error of Az < 0.003(1 + z), volume of V' = lOGpc'^/h^ 
and number density of fi = 10~^/i''/Afpc^ LRG galaxies. 
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